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Abstract 


A  new  numerical  method,  “Wave  Confinement”  (WC),  is  developed  to  efficiently  solve  the  linear 
wave  equation.  This  is  similar  to  the  originally  developed  “Vorticity  Confinement”  method  for 
fluid  mechanics  problems.  It  involves  modification  of  the  discrete  wave  equation  by  adding  an 
extra  nonlinear  term  that  can  accurately  propagate  the  pulses  for  long  distances  without  numerical 
dispersion/d ifFus ion.  These  pulses  are  propagated  as  stable  codimension-one  surfaces  and  do  not 
suffer  phase  shift  or  amplitude  exchange  in  spite  of  nonlinearity.  The  pulses  remain  thin  unlike 
conventional  higher  order  numerical  schemes,  which  only  converge  as  N  (number  of  grid  cells 
across  the  pulse)  becomes  large.  The  additional  term  does  not  interfere  with  conservation  of  the 
important  integral  quantities  such  as  total  amplitude,  centroid.  Also,  properties  like  varying  index 
of  refraction,  diffi^tion,  multiple  reflections  are  included  and  tested. 

The  generated  short  pulses  can  be  best  described  as  solitary  waves,  which  can  recover  the  shape 
after  a  collision  due  to  nondestructive  interaction  between  the  pulses.  Within  the  pulse,  the 
dissipative  effects  due  to  the  numerical  errors  are  balanced  with  those  of  nonlinearity  and  the 
pulse  will  its  their  original  form  and  speed  even  after  many  collisions.  The  pulse  is  also  used  as  a 
carrier  wave  to  propagate  other  properties  such  as  direction.  Wave  equation  solutions  in  the  high 
frequency  approximation  can  be  generated  using  the  carrier  wave  approach.  WC,  together  with 
Keller's  Approximation  is  then  used  to  capture  diffraction  effects  from  a  straight  edge. 

Scattering  over  complex  bodies  can  be  modeled  with  no  use  of  complicated  adaptive  grid 
generation  schemes  around  the  bodies.  The  confinement  term  smoothens  the  boundary  and 
prevents  stair  casing  effects  but  the  boundary  remains  thin.  Validation  studies  have  been 
performed  for  a  number  of  real  flow  models  and  compared  to  the  exact  solutions.  It  is  observed 
that  the  solutions  match  quite  well  with  the  exact  solution. 

1.  INTRODUCTION 

The  main  objective  of  this  dissertation  is  to  develop  a  new  and  efficient  algorithm  to  accurately 
solve  the  linear  wave  equation.  The  idea  comes  from  an  already  existing  method  “Vorticity 
Confinement”  (VC),  which  was  developed  for  a  vast  range  of  fluid  dynamics  problems  [ij. 
Hence,  the  new  method  “Wave  Confinement”  (WC)  is  named  after  the  existing  approximation 
(VC).  The  linear  wave  equation  can  describe  different  classes  of  wave  phenomena.  Such  classes 
include  acoustics,  electro-magnetics,  microwave  theory,  etc.  This  dissertation  focuses  on  scalar 
wave  propagation  over  long  distances.  Effects  to  be  accommodated  include  variation  of  the  index 
of  refraction,  multiple  scattering  from  complex  surfaces,  and  some  cases  of  diffraction. 


2.  APPROACH 
2.1  Introduction 

Let  us  consider  a  scalar  field,  which  satisfies  the  one  sided  wave  equation  (advection 
equation)  advecting  with  constant  speed  c: 


2.1 


— i-  +  C  — 
dt  dx 


=  0 


Evidently,  if  the  above  equation  is  numerically  solved  using  a  conventional,  discrete  finite 
difference  scheme,  it  tends  to  develop  discretization  errors.  The  conventional  Eulerian  schemes 
[8]  suffer  dissipation  no  matter  what  the  order  is.  Many  higher  order  schemes  have  been  proposed 
to  reduce  the  numerical  diffusion/dispersion,  but  they  only  reduce  the  error  if  the  numb^  of  grid 
points  across  the  pulse,  N  is  relatively  large.  As  the  present  work  involves  treating  thin  pulses 
2-3  grid  cells  in  size,  the  method  is  not  necessarily  going  to  be  more  accurate  with  increasing 
order.  To  keep  the  pulse  confined  on  the  discrete  domain,  numerical  dissipation/dispersion  must 
be  counter-acted  by  adding  a  new  term  to  the  advection  equation  that  will  not  interfere  with  the 
essential  properties  of  the  pulse.  This  term  is  called  the  “Confinement”  term  and  hence  the 
method  is  named  “Wave  Confinement”  (WC).  Ideas  of  adding  extra  terms  to  prevent  numerical 
dissipation  were  previously  proposed  by  Harten  [9]  to  capture  ID  contact  discontinuities  in 
compressible  flows.  The  discretized  form  can  only  conserve  a  limited  number  of  physical 
quantities  of  the  pulse  corresponding  to  the  limited  number  of  grid  nodes  across  the  pulse.  Here, 
these  are  taken  to  include  total  amplitude  at  each  grid  point  and  the  centroid  speed.  The  modified 
form  of  equation  (2. 1 )  with  the  confinement  term  is  then 


5/  dx 


2.2 


where,  E  =  — — .  The  idea  is  that  the  second  order  form  of  F  acts  as  a  “pulse  shaping”  term  and, 
dx 

as  long  as  F  — >  0  as  jx|  oo ,  the  quantities  of  interest  are  conserved.  A  simple  form  of  F  can  be 


taken  as 


2.3 


where  //' ,  .  Here,  constants  and  O  is  the  harmonic  mean  given  as 


where  h  is  the  distance  between  two  points,  which  tends  to  zero  in  the  continuous  limit.  The 
confinement  term  is  chosen  such  that  it  acts  as  an  “expansion”  or  positive  diffusion  for  shorter 
wavelengths  and  a  “contraction”  or  negative  diffusion  for  longer  wavelengths.  Also,  the  method 
must  have  a  nonlinear  term,  otherwise  some  modes  will  decouple  and  diverge.  The  nonlinear  term 
prevents  the  amplitude  to  escape  to  long  wavelengths  and  can  keep  the  pulse  confined  to  a  few 
grid  cells. 


2.2  WC  as  a  Nonlinear  Partial  Differentia]  Equation 


To  better  understand  the  properties  of  the  WC  method,  the  harmonic  mean,  O,  given  in  equation 
(2.4)  is  approximated  as  a  partial  differential  equation.  The  Taylor  expansions  for  + A)  and 
^(jf  -  h)  to  second  order  are 


f>(jc  + ;»)  =  ^{x)  +  hl)\x)  +  ^  f{x) , 


2.5 


where  =  =  By  substituting  the  Taylor  expansions  into  equation  (2.4),  and 

neglecting  the  higher  order  terms,  the  nonlinear  term,  O ,  is  then  approximated  as 


2.6 


where  A  — >  0  in  continuous  limit  and  terms  of  the  order  can  be  neglected.  Using  equation 
(2.6),  equation  (2.2)  becomes 


2.7 


When  =  the  above  equation  fails  but,  it  is  important  to  note  that  the  discretized  form  of 
equation  (2.7)  used  for  computations  is  given  in  equation  (2.17).  The  harmonic  mean  in  equation 
(2.2 1 )  is  still  finite  for  ^  =  0  . 

Equation  (2.7)  describes  the  conservation  of  the  following  integral  quantities: 
total  amplitude. 


2.8 


and  centroid  speed, 


2.9 


dt  Af 


where  x  is  the  centroid,  given  by 


2.10 


It  IS  important  to  note  that  all  the  terms  in  equation  (2.7)  are  homogenous  of  degree  one  unlike 
many  non<linear  equations,  which  use  non-homogenous  terms  for  the  nonlinear  term  [10].  Tlius, 
the  confinement  does  not  depend  on  the  scale  of  the  quantity  to  be  convected  (which  is,  of  course, 
a  property  of  the  original  linear  equation  (2.1)).  Defining  the  confinement  term  in  the  equation 
(2.7)  as  the  sum  of  three  terms  gives, 


£  =  +  £,  +  £2 , 


2.11 


second  order  term  (  Eq)  in  equation  (2.1 1)  is  different  from  typical  nonlinear  pde’s  studied,  such 
as  KdV,  that  harbor  solitary  waves:  In  these,  the  linear  second  order  term  is  the  “expansion”  term, 
and  the  “contraction”,  or  “steepener”  term  is  the  nonlinear  Burgers-like  convection:  (5^^^  /2). 
In  WC,  the  linear  second  order  term,  Eq,  contracts  the  pulse  and  the  nonlinear  term,  £3, 


prevents  ^  from  changing  sign  and  transfers  amplitude  from  long  wavelength  to  short,  and  the 


fourth  order  term,  £, ,  acts  as  a  diffusion  for  short  wavelengths  and  prevents  the  solution  from 
diverging.  The  amplitude  imparted  to  the  field  by  the  WC  terms  remains  confined  and  propagated 
as  stable  wave  packets  without  dispersion/diffijsion.  The  appearance  of  ^  in  the  denominator  of 
equation  (2.7)  makes  £2  diverge  as  »0.  This  prevents  ^  from  changing  sign.  Since  the  total 
amplitude,  ,  is  conserved,  the  integral  of  ^  over  any  finite  region  cannot  then  diverge.  In  the 
discretized  version  defined  in  Section  2.3,  none  of  the  grid  values  can  diverge  in  appropriate 
ranges  of  a  and  A.  This  ensures  realizability  if  ^  is  a  physical  quantity.  Smolarkiewicz  [11] 
also  has  rearranged  the  discretized  convection  equation  so  that  there  is  such  a  term  in  the 
denominator  for  this  reason. 

By  rearranging  the  terms,  equation  (2.7)  may  be  expressed  as 


2.12 


where  In  the  convecting  frame,  ^  =  x - cf ,  the  above  pde  simply  reduces  to  the  heat 

equation. 


d,^  =  d]F 

where  F  =  — ^  (—  aif/  +  At  convergence,  F  — >  0 ,  and  therefore  reduces  to 


2.13 


-ay/  +  Xd\y/  =  0. 


2.14 


It  can  be  seen  that  the  above  equation  is  similar  to  the  simplest  form  of  Sturm-Liouville 
eigenvalue  equation  with  P  -  yJcc/A ,  The  solution  to  the  equation  (2.14)  will  then  be  of  the 
form, 

^  =  C,(e^'  2.15 


Then  the  solution  for  ^  is 


2.16 


where  C, ,  Cj  and  A  are  arbitrary  constants. 

2.3  Discretization 

Using  first  order  upwind  in  time  and  second  order  centered  in  space,  the  discretized  form  of 
equation  (2.2)  is  written  as 


2.17 


rAt 

where  SjFj  =  ^  ^  the  time  step  and  h  is  the  grid  cell  size.  Even 

on  the  discrete  domain,  the  confinement  term,  F,  conserves  the  essential  quantities  such  as  total 
amplitude, 


Ar  =  const. 


2.18 


and  centroid  speed. 


Ar 

where  J  is  the  total  number  of  grid  points.  The  confinement  term  F  is  defined  as 


2.19 


2.20 


Z/'A/  ,  , 

where  //  =  — — ,  S  =  — — ,  /z  is  a  diffusion  coefficient  and  S  is  the  confinement  coefficient. 

h  h 

E*  along  with  /z' ,  control  the  width  and  rate  of  convergence  of  the  pulse,  and  is  the  nonlinear 
harmonic  mean  of  N  (including  N-l  neighboring  nodes)  grid  points.  For  1 D,  N  =3,  and, 


2.21 


Equation  (2.1 7)  would  simply  reduce  to  the  Euler  explicit  form,  which  is  unconditionally  unstable 
when  s  =  0  and  /i  =  0 .  The  two  (positive)  parameters,  s  and  /i ,  are  determined  by  the 
requirements  for  convergence,  such  that  the  small  features  relax  to  their  solitary  wave  shape  in  a 
small  number  of  time  steps  and  can  have  an  effective  support  (to  have  a  significant  magnitude)  of 
a  small  number  of  grid  cells.  In  the  semi-discrete  limit,  the  scalar  field  ^  relaxes  to  a  definite 
shape  as  derived  in  Section  2.3.1.  In  the  moving  frame  (v'  =  0),  equation  (2.17)  reduces  to  a 
discrete  heat  equation, 


r;'=rj+5]F';  2.22 

that  simulates  the  stationary  pulse.  At  convergence,  the  positive  and  negative  artificial  dissipation 
terms  are  balanced  with  the  nonlinearity  to  produce  stable  pulses  that  can  stay  confined  and 
structurally  stable  to  numerical  perturbations  (caused,  for  example,  by  the  discrete  convection 
term  or  another  wave  passing  through). 

Zi./  Convergence 

At  convergence,  SjF  F  ->0  and  (//fJ-fO)  — ►  0 .  IfN  =  3  in  the  harmonic  mean,  the 
equation  for  F  at  a  point  is 


2.23 


Solutions  to  equation  (2.23)  that  satisfy  the  boundary  conditions:  F  0,  |yj  oo  are  of  the 
form. 


4>^A%tch[YU-jo)]> 


2.24 


where  is  the  approximate  position  of  the  centroid.  Here,  is  a  pulse  width  coefficient  and  is 
determined  by  substituting  equation  (2.24)  into  (2.23): 


coshy  = 


2.25 


Furthermore,  the  above  equation  can  be  reduced  to  a  quadratic  form  in  e^  ^ 


for  which  the  solution  is 


6^  = 


b±  -4 

2 


2.27 


or 


r 


/ 


2.28 


where  b  = 


\ 

1  .  The  solution  is  real  if  6^  -  4  >  0  and  this  requires  that  s//i^].  If  f  =  0 , 


F  acts  as  a  positive  diffusion  for  all  wavelengths  and  the  solution  decays.  When  //  =  0,  the 
solution  becomes  unstable  even  if  f  =  0.  So,  a  minimal  amount  of  diffusion  is  essential  for  the 
stability  purpose.  If  ^  ,  F  acts  as  negative  diffusion  for  longer  wavelengths  so  that  the 
features  remain  confined  and  stable  to  perturbations  against  spreading. 


2.4  Results 

For  the  evaluation  of  the  scheme  and  to  demonstrate  the  efficient  properties  of  the  WC  method,  a 
number  of  tests  are  performed  and  presented.  To  show  convergence,  a  single  point  pulse  of 
arbitrary  amplitude  chosen  to  be  2, 


=  2.29 

lO.  J*Jo\ 

with  parameters  »^  =  0,  f  =  0.3,  //  =  0.2  and  =  128  is  simulated  on  a  grid  of  256  points 
using  pa*iodic  boundary  conditions.  The  discretized  equation  used  for  this  computation  is 

+  2.30 

When  the  solution  reaches  convergence,  the  pulse  is  relaxed  to  a  hyperbolic  secant  given  by 
equation  (2.24).  In  Figure  2.1,  the  solid  line  is  the  algebraically  calculated  hyperbolic  secant  with 
width  given  by  equation  (2.28),  and  the  points  are  computed  using  equation  (2.30).  For  better 
comparison,  the  normalized  function,  A  is  plotted  against  the  hyperbolic  secant.  It  can  be 
clearly  seen  that  the  solution  matches  quite  well  with  the  analytically  calculated  hyperbolic 
secant.  When  discretized  form  of  equation  (2.1)  is  solved  using  a  higher  order  method,  the 
solution  quickly  spreads  to  a  large  number  of  grid  points.  It  can  be  seen  in  Figure  2.2  that  the 


solution  from  a  higher  order  method  spreads  by  large  amount  af\er  only  100  times  steps,  but  the 
confinement  method  keeps  the  pulse  compact. 

3.  FORMULATION 

As  in  convection,  a  nonlinear  term,  E  is  added  to  the  one  dimensional  full  wave  equation.  The 
wave  equation  with  confinement  terms,  that  control  the  shape  of  the  pulse,  is 


d^(fi  =  c^d]^  +  E  3.1 

where  E  =  d,d\F .  A  simple  discretization  for  the  above  pde  will  be 

f/'  =  Ifj  -  fj-'  +  +  S;S]F  3.2 

cA/ 

where  ,  S^fj  =  -  2fj  +  ,  v  =  — —  ,  A/  is  the  time  step  and  h  is  the 

h 

grid  cell  size,  and  F  is  defined  as  in  the  advection  case. 


3.3 


where  fi  = 


/t’A/ 


E  = 


ff'A/ 


solution  to  the  equation  (3.2)  is 


.  At  convergence  (when  the  propagating  pulse  shape  converges),  the 


(p^A^  sec  h{/{J  -  +»«))+ j  sec  hiyiJ  -j^-vn )) 


3.4 


where  and  are  arbitrary  constants,  is  the  initial  position  of  the  centroid  and  the  pulse 
width  coefficient,  y ,  is  defined  as  in  equation  (2.28).  Equation  (3.4)  describes  two  pulses  moving 
forward  and  backward.  It  is  shown  in  Section  2  that  the  addition  of  WC  terms  in  the  form  of 
second  derivatives  of  a  function  that  has  short  range  do  not  change  the  propagation  speed  or  the 
total  amplitude.  The  same  is  true  for  the  wave  equation,  if  an  additional  time  derivative  is  applied. 
The  main  constraint  on  the  confinement  term,  F,  as  in  advection  is  that  it  forces  an  initial  isolated, 
propagating  short  range  pulse  with  single  maxima  to  remain  short  range  and  also  not  develop  any 
additional  maxima.  An  important  property  of  the  wave  equation  is  that  it  is  linear.  Therefore, 
properties  of  the  computed  solution  should  not  depend  on  the  amplitude.  The  behavior  of  solution 
computed  using  equation  (3.1)  is  shown  in  Figure  3.1.  When  //-O  and  f  =  0,  the  solution 
becomes  unstable.  A  little  amount  of  diffusion  (//  =  0.2)  will  prevent  the  unstable  behavior  but 
the  solution  suffers  large  amount  of  spreading.  Now,  when  the  confinement  term  (£*=0.3)  is 
added,  the  pulse  becomes  stable  and  at  the  same  time,  remains  thin.  The  dissipative/dispersive 
effects  are  balanced  with  those  of  the  nonlinearity  to  produce  stable  localized  structures  [13]. 


3.2  Results 


In  this  section,  numerical  examples  are  provided  to  confirm  that  no  discretization  errors  are 
developed  in  the  quantities  of  interest  and  that  there  is  no  significant  spreading  of  the  pulses 
during  propagation.  The  formulation  given  by  equation  (3.4)  is  used  for  the  numerical 
experiments.  While  analyzing  the  formal  accuracy  of  the  present  scheme,  a  theoretical  error 
estimate  is  beyond  the  scope,  because  the  equation  is  nonlinear.  However,  evaluation  of 
numerical  error  is  done  experimentally,  and  demonstrated  that  the  scheme  possesses  the 
convergence  that  was  discussed  in  the  previous  chapter. 


3,Z1  Long-time  propagation 


To  begin  with,  the  propagation  of  pulses  for  long  distances  is  observed.  This  is  the  idealization  of 
waves  without  dissipation.  Convergence  analysis  is  illustrated  by  examining  the  behavior  of  a 
delta  function  as  an  initial  pulse.  Let  the  initial  pulse  be 


J=jo 

j^Jo 


at  n  =  0  and  n  =  I 
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where  =  1 28  and  the  computation  is  done  on  256  grid  points.  The  pulse  from  WC  is  compared 
to  the  solution  from  a  conventional  higher  order  method  and  is  shown  in  Figure  3.2.  The  two 
pulses  in  the  plot  are  the  forward  and  backward  propagating  scalar  functions.  It  is  obvious  that 
WC  preserves  the  thickness  and  total  amplitude  of  the  pulse  with  no  significant  numerical 
dissipation.  Higher  order  methods  do  not  help  because  increasing  the  order  of  the  numerical 
scheme  is  only  efficient  when  there  are  a  large  number  of  points  across  the  propagating  function. 

3.2.2  Wave  Interaction 


Before  moving  to  higher  dimensions,  the  pulse  interaction  phenomenon  is  discussed  to 
understand  the  performance  of  WC  during  multiple  collisions  of  pulses.  To  represent  interacting 
wave  equation  solutions,  when  pulses  pass  through  each  other,  there  must  be  no  amplitude 
exchange  or  phase  shift.  Otherwise,  the  actual  wave  equation  being  studied  could  not  be 
accurately  simulated,  since  it  is  linear.  However,  a  nonlinear  term  is  required  in  the  WC  equation 
in  order  to  create  a  solitary  wave  representation  of  the  pulse,  which  will  be  non-diffusing  when 
discretized.  Interaction  between  soliton-like  pulses  is  an  important  property,  which  has  been 
studied  for  a  long  time  using  nonlinear  equations  [14].  It  is  more  appropriate  to  only  re-introduce 
the  nonlinear  pde  for  WC,  which  was  derived  in  Section  2,  to  establish  a  relation  to  other  existing 
nonlinear  pde’s.  The  nonlinear  pde  for  WC  is 


d,ij>=d\ 


2\ 
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The  above  equation  has  a  similarity  to  the  Cahn-Hilliard  (CH)  equation,  which  describes  the 
process  of  phase  separation,  by  which  the  two  components  of  a  binary  fluid  spontaneously 
separate  and  form  domains  pure  in  each  component  A  simple  form  of  the  CH  equation  is 

(-  «,<*  -  a^d\4)  +  a,4>^) 
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where  a, ,  a, ,  a,  arc  arbitrary,  positive  constants.  The  commonality  includes  the  same  linear 
terms  and  the  appearance  of  a  Laplacian  in  front  of  the  non-linear  term.  This  commonality  is  not 
surprising  since  the  CH  equation  was  one  of  the  leading  models  proposed  for  a  phenomenological 
description  for  a  fluid  interface  and  phase  separation  (decomposition  into  pure  phases)  [15].  The 
CH  equation  describes  the  evolution  of  concentrated  fields  (like  thin  pulses  evolved  from  WC). 
Though  the  equations  are  similar,  the  resulting  non-linear  term  is  not  exactly  the  same  as  in  the 
commonly  used  form  of  the  CH  equation. 

An  important  point  is  that  other  nonlinear  pde’s  like  Kdv,  which  can  successfully  propagate 
localized  structures,  suffer  phase  shift  after  collision  [16].  This  is  not  present  in  WC.  For  the  new 
formulation  used,  this  interaction  effect  vanishes.  This  vanishing  is  due  to  the  fact  that  both  the 
Laplacian  and  the  time  derivative  operator  operate  on  the  nonlinear  term. 

The  preservation  of  centroid  speed  during  interaction  is  shown  in  Figure  3.3  for  an  initial 
condition  given  by  equation  (3.5).  The  centroid  positions  of  forward  and  backward  propagating 
pulses  are  plotted.  Periodic  boundary  conditions  are  used  for  this  simulation  over  a  grid  of  256 
grid  points.  It  can  be  noticed  that  even  after  many  collisions,  no  jump  in  the  centroid  position  is 
noticed  after  interaction.  The  importance  of  preserving  the  right  speed  plays  a  dominant  role  in 
generating  constant  phase  surfaces,  which  will  be  discussed  in  Section  5.  The  same  is  true  for 
preservation  of  total  amplitude.  In  Figure  3.4,  it  can  be  seen  that  two  pulses  of  different 
amplitudes  are  effectively  transparent  (after  a  short  relaxation  time)  to  one  another  and  do  not 
lose  their  amplitudes  during  collisions.  For  this  computation,  the  initial  conditions  used  are. 


J  ~  Jo\ 

fj  1.  7=7*02  ^  3.9 

0,  7  ^  7oi  7  ^  7o2 

The  two  initial  pulses  are  first  relaxed  to  hyperbolic  secants  of  respective  amplitudes  using  heat 
equation  (c->(?).  So,  the  two  initial  conditions  of  the  two  pulses  are  then 


(a,  sec//[^(y  -70,)]  +  ^  sec//[y'(y-  yoi)],  n  =  0 

\^sech[r{j- Jo,-v)]+  A^sechlriJ- Jai+y)],  n  =  \ 


where  /Ij  =  =  128  and  Jq2  =  140 .  The  pulses  are  then  propagated  towards  each  other  to 

observe  the  interaction  between  the  pulses. 

33  WC  for  higher  dimensions 
3JJ  2D 


In  higher  dimensions,  WC  behaves  the  same  way  as  in  ID.  The  essential  physics  is  accurately 
preserved  as  analyzed  in  Section  2.  The  discretized  equation  in  2D  is 


3.11 


where  /  is  the  grid  index  in  x-dircction,  y  is  the  grid  index  in  y-direction  and  is  the  2D 
extension  of  the  discrete  Laplacian  used  in  ID  .  For  convenience  the  aspect  ratio  of  the  grid  is  1. 
However,  behavior  of  WC  with  different  aspect  ratios  is  studied  in  Section  4.  The  nonlinear 
factor,  ,  in  the  confinement  term  F  is 

•J"  N 


where  the  number  of  terms  in  the  sum,  =  5.  Here,  it  is  assumed  that  >  0 .  Negative  values 
can  also  be  accommodated  with  a  small  extension.  Both  /i  and  €  are  positive.  Assuming 
convergence  as  (for  K  =  0  ), 


=  0  3.13 

for  which  the  solution  at  convergence  is 

^,j=AsccHy(r,j-rg)]  3.14 

where 

r,j  =XfCx>sS  +  yjS\nS  3.15 

A,  Zq,  9  are  constants,  jc,  =/A,  =  Jh^  h  is  the  grid  cell  size,  and  y  is  the  inverse  pulse 

width.  The  solution  converges  to  a  straight  pulse  (in  2-D)  concentrated  about  a  line  at  angle,  9 
and  at  any  position.  By  substituting  the  solution  in  equation  (3.13),  it  is  easy  to  see  that  / 
satisfies 

—  =  [\'\-2ch(yhoos9)  +  2ch(}hsin9)]/5.  3  16 

M 

The  solution  has  translational  invariance  and  y  depends  (to  a  small  extent)  on  orientation,  9, 
However,  the  centroid  speed  and  amplitude  in  the  normal  direction  is  independent  of  ^  as  can  be 
seen  from  equation  (3.16).  To  demonstrate  long  distance  propagation  with  no  increasing 
numerical  errors,  a  diverging  circular  wave  of  speed,  v  =  0.23  is  simulated  on  a  (128)^  cell  grid 
with  periodic  boundary  conditions.  Confinement  values  used  were  /;  =  0.2 ,  £■  =  0.3 .  Tbe  initial 
conditions  for  this  computation  are 
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where  r  is  the  radius  of  the  wave.  The  actual  wave  equation  exhibits  a  ‘"tail”  behind  a  pulse  in  2- 
D,  which  can  be  seen  to  be  suppressed  by  the  Confinement,  and,  effectively,  only  the  steep  pulse 
front  is  accurately  computed.  In  Figure  3.5,  the  circular  wave  for  different  time  steps  is  presented. 
It  can  be  seen  that  the  wave  persists  on  the  grid  for  veiy  long  time.  Also,  as  in  I  -D.  there  is  no 
discemable  interaction  between  intersecting  waves.  The  waves  retain  their  form  and  orientation  in 
spite  of  multiple  head-on  collisions. 

33.2  3D 


In  3D,  the  discretized  form  of  the  wave  equation  with  the  confinement  term  is 


K*j\  •  Kj.*  -C + ^  ) 
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where  is  the  3D  extension  of  the  discrete  Laplacian  used  in  ID,  and  the  nonlinear  harmonic 
mean  at  a  point,  *  ,  is 


<t>"  .  = 


4-1  *]  *\ 


cr— 1 Ijr—t _ 

27 
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Long  distance  propagation  in  3D  is  simulated  for  an  expanding,  initially  spherical  wave  and  is 
shown  in  Figure  3.6.  The  computation  is  done  on  a  coarse,  (64)^  cell  grid  with  periodic 
boundary  conditions.  The  initial  diameter  for  this  computation  is  32  grid  cells.  Confinement 
values  used  were  ^  =  0.2  and  e  =  0.3 .  As  in  the  2-D  case,  the  wave  remains  completely 
confined.  Unlike  in  2D,  the  actual  solution  does  not  exhibit  a  tail  behind  the  front  and  the  solution 
corresponds  to  the  full  physical  solution  rather  than  just  the  wavefront.  The  solution  at  a  given 
point  P(x,  y,  z)y  depends  only  on  the  information  on  the  sphere  and  not  on  the  interior  of  the 
sphere.  Thus  the  interior  of  the  sphere  is  a  lacuna  for  the  solution.  It  is  true  for  odd  numbers  of 
space  dimension. 


4  WAVE  FRONT  PROPAGATION 
4.1  Focusing 

Focusing  is  the  phenomenon  of  the  waves  converging  on  a  point  or  line.  This  is  important  in 
studying  caustic  regions.  A  common  situation  where  caustics  are  visible  is  when  light  shines  on  a 
drinking  glass.  Focusing  phenomenon  has  been  extensively  applied  in  a  variety  of  engineering, 
optics,  computer  graphics  and  medical  applications.  In  caustic  regions,  many  waves  intersect. 
Fixed  grid  Eulerian  methods  have  limits  on  the  effective  resolution  due  to  the  limited  number  of 
grid  points  in  a  caustic  unlike  Ray  Tracing  schemes  [4]  [17].  To  overcome  this  problem,  complex 


schemes  have  been  considered  such  as  Osher’s  higher  dimensional  Liouville  equation 
representation  [18]  and  localization  of  caustics  [19].  However,  these  seem  to  be  complex  and 
computationally  costly. 

It  was  first  thought  that,  since  during  focusing,  the  number  of  grid  points  on  the  wave  decreases, 
WC  is  not  going  to  help  to  retain  the  information  at  focusing  regions.  It  was  later  shown  that  WC 
accurately  computes  waves  through  caustic  regions  and  is  stable  to  discretization  effects  produced 
on  the  grid.  Since  the  interest  is  in  the  long  distance  propagation,  the  detailed  resolution  at  the 
focusing  itself  was  not  an  issue  and  that  focusing  in  interniKdiate  regions  docs  not  interfere  with 
the  long  distance  wave  propagation.  In  Figure  4.1,  amplitude  contours  are  compared  to 
Lagrangian  markers  (exact  solution  in  high  frequency  approximation).  It  can  be  seen  that  the 
wave  contours  match  quite  well  with  markers.  Also,  the  basic  information  defined  by  the  initial 
conditions  is  not  lost,  even  though  only  the  simple  Eulerian  algorithm  was  used,  with  no 
additional  logic.  Recent  developments  in  level  set  methods  can  capture  focusing  regions  [21]. 
However,  these  are  computationally  complex  unlike  WC. 


4.2  Reflection  from  Complex  Boundaries 

The  property  of  waves  to  reflect  from  boundaries  is  captured  accurately  using  WC  without  having 
to  use  adaptive  grids  around  the  boundaries.  This  reflection  of  waves  is  responsible  for  echoes, 
interference,  etc.  The  energy  and  momentum  of  the  waves  is  reflected  back  with  or  without  1 80 
degrees  phase  change  depending  on  the  properties  of  the  boundary.  The  form  of  a  reflected  wave 
front  is  determined  by  those  of  the  incident  wave  front  and  the  surface.  Boundaries  many  contain 
many  irregularities  and  for  short  pulses  (high  frequency  content),  these  irregularities  can  be  larger 
than  the  wavelengths.  Higher  order  discretization  may  be  needed  for  conventional  treatments  of 
long  distance  propagation  to  reduce  numerical  diffusion  in  the  propagation.  Since  WC  does  not 
have  this  diffusion  problem,  complex  configurations  can  be  accurately  captured  using  very 
efficient,  lower  order  discretizations  with  no  loss  of  accuracy.  A  very  effective  method  for 
treating  reflections  can  then  be  implemented  that  docs  not  require  complex  surface  fitted  grids, 
but  allows  the  surface  to  be  simply  ‘immersed’  in  a  uniform  Cartesian  grid  as  shown  in  Figure 
4.2.  This  method  employs  a  “level  set”  representation  of  the  surface  and  can  easily  accommodate 
very  complex  topography  with  little  computational  effort.  For  example,  reflection  from  a 
boundary  of  function,  fiUj)  =  0  ,  is  shown,  where, 

=  0,  /  ^  0,  at  all « ,  4.1 

The  effects  of  boundaries  are  smooth  and  no  stair  casing  is  observed.  Reflection  from  complex 
terrains  with  many  irregularities  is  important  to  capture  interference  effects  in  radio  wave 
propagation.  The  irregularities  can  be  very  much  larger  than  the  wavelength  of  the  short  of  waves 
that  are  propagated. 

4.2.1  Scattering  from  Cylinder 

Scattering  is  a  physical  process  in  which  the  direction  of  a  wave  deviates  from  its  normal 
trajectory  due  to  non-uniformities  in  the  medium  or  interaction  with  boundaries.  A  wave  packet 
with  a  bundle  of  rays  parallel  to  each  other  will  reflect  in  different  angles  if  they  encounter  a  non- 


uniform  obstacle  (diffose  reflection  or  scattering).  If  such  a  bundle  encounters  a  smooth  surface 
such  as  mirror,  the  rays  remain  as  bundle  upon  leaving  the  surface  (specular  reflection). 

Scattering  from  the  cylinder  is  diffuse  reflection  and  is  simulated  accurately  using  WC.  Consider 
the  oblique  incidence  of  a  plane  wave  on  a  perfectly  conducting  infinite  cylinder.  The  scalar  field, 

(f> ,  is  set  to  zero  inside  the  cylinder  of  radius  32  grid  cells.  The  reflected  surfoce  is  manipulated 
to  re-direct  or  re-shape  the  plane  wave  front  into  a  cylindrical  scattered  field.  If  the  plane  wave  is 
polarized,  the  polarity  is  reversed  when  it  encounters  perfectly  conducting  obstacles  or  hard 
boundaries  (This  is  not  demonstrated  here).  The  confinement  term,  F  is  diffusive  in  the  tangential 
direction  and  contracting  in  the  normal  direction.  This  prevents  stair-casing  effects  and  generates 
a  smooth  scattered  field,  which  can  be  seen  in  Figure  4.3.  The  scattered  field  does  not  quickly 
diffuse  on  the  grid  and  can  be  detected  near  the  receiver,  which  is  located  far  away  from  the 

target.  This  is  important  in  applications  such  as  target  detection.  The  scalar  quantity,  ^ ,  is  very 
sensitive  and  can  be  used  to  detect  targets  smaller  than  the  grid  cell  size.  When  the  field  of 
computation  is  very  large,  the  number  of  grid  points  across  the  target  will  be  relatively  small.  The 
direct  scattering  problem  involves  determining  the  scattered  field,  which  depends  on  the 
characteristics  of  the  scatterer  studied  in  this  section  for  different  cases.  However,  a  lot  of  work  is 
yet  to  be  done  in  inverse  scattering,  which  is  considered  for  future  research. 

Before  proceeding  to  the  study  of  complex  configurations  of  obstacles,  a  plane  wave  reflection 
from  a  sine  wave  boundary  is  computed  and  shown.  The  same  can  be  effectively  captured  using 
WC  and  is  shown  in  Figure  4.4.  It  can  be  seen  that  WC  itself  can  capture  the  effects  of  the 
boundary. 


4.3  Varying  Index  of  Refraction 

When  a  wave  interacts  with  the  medium  of  varying  density,  it  refracts  (scatters).  Refraction  can 
be  described  by  Snell’s  law,  which  states  that  the  angle  of  incidence  is  related  to  the  angle  of 
refraction. 


sin  K 

=  — =  —  4.2 

Sin  ^2  ^2 

where  n, ,  are  the  indices  of  refraction  of  the  two  layers  of  medium,  6^ ,  $2  are  the  angles  of 
incidence  and  refraction  respectively  and  v, ,  1/2  are  the  speeds  of  the  pulse  in  the  respective 
layers.  One  example  of  refraction  in  the  atmosphere  is  mirage. 

In  acoustics,  refraction  is  the  bending  or  curving  of  a  sound  ray  that  results  when  the  ray  passes 
through  a  sound  speed  gradient  from  a  region  of  one  sound  speed  to  a  region  of  a  different  speed. 
The  amount  of  ray  bending  is  dependent  upon  the  amount  of  difference  between  sound  speeds. 
For  water,  sound  speed  depends  on  the  variation  in  temperature,  salinity,  and  pressure  of  the 
water.  Similar  acoustics  effects  are  also  found  in  the  Earth's  atmosphere.  The  phenomenon  of 
refraction  of  sound  in  the  atmosphere  has  been  known  for  centuries;  however,  beginning  in  the 
early  1970s,  widespread  analysis  of  this  effect  came  into  existence  through  the  designing  of  urban 
highways  and  noise  barriers  to  address  the  meteorological  effects  of  bending  of  sound  rays  in  the 


lower  atmosphere.  Refniction  of  radio  waves  in  evaporated  ducts  (formed  due  to  variations  in 
humidity  near  sea)  and  ionospheric  ducts  (formed  due  to  variation  in  solar  radiation)  is  important 
in  radio  wave  propagation. 


4,3,1  Profile  of  varying  index  of  refraction:  Converging  waves 
Consider  a  profile  for  varying  speed, 

v{j)  =  const  *cxp(aU-Joy),  4.3 


which  has  maximum  speed  (i^ )  at  J  -  Jq-  An  isolated  plane  wave  propagated  across  the  medium 

with  the  profile  given  in  equation  (4.3)  is  computed.  The  discretized  equation  used  for  the  wave 
propagation  is, 


c = Kj  -  c/  *  ^yjKj + syjKj + + s]  ) . 
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When  the  index  of  refraction  is  non-uniform,  the  waves  tend  to  bend  towards  regions  of  higher 
index.  It  can  be  seen  in  Figure  4.8  that  the  computed  pulses  travel  without  diffusion  and 
numerical  errors  in  speed  were  insignificant  to  plottable  accuracy.  The  computed  pulse  is 
compared  to  Langrangian  markers  (the  blobs  in  Figure  4.5)  for  validation,  which  are  calculated 
as, 


4.5 


where  and  are  the  components  of  the  propagating  direction.  The  direction  is  updated 
according  to  the  gradient  in  index  of  refraction  on  that  marker  as. 


dv 

dy 


and  the  normalized  components  of  the  direction  are. 


^x* 


''x^new 


+  v' 


4.6 


4.7 


When  the  index  of  refraction  is  randomly  vaiying  and  does  not  have  a  smooth  profile,  curve 
fitting  methods  must  be  used  to  fit  the  random  variations.  This  becomes  very  complex  when  the 
index  of  variation  is  a  f\inction  of  time.  But,  for  a  smooth  profile,  the  computation  of  markers  is 
exact  and  is  used  for  validation. 

As  shown  in  Figure  4.5,  the  pulse  does  not  lose  information  in  spite  of  the  limited  number  of  grid 
points  across  the  focusing  regions.  However,  Lagrangian  markers  suffer  inadequacy  in  number  of 
grid  points  and  fail  to  maintain  continuity. 

4.3,2  Duct  like  profile  used  in  EM  propagation 

Electromagnetic  waves  propagating  in  the  atmosphere  encounter  regions  of  varying  density, 
which  cause  bending  of  the  waves.  One  example  is  the  evaporation  duct,  which  is  associated  with 
a  sharp  drop  in  moisture  immediately  above  water  bodies.  The  drop  in  moisture  decreases  the 
density  and  in  turn  increases  the  speed. 

A  predetermined  profile  of  varying  index  of  refraction,  which  imitates  such  a  duct  is  treated: 

Ky)  =  0.23- 0.0023  ♦|7-7o|  4.10 

where  Jq  =  30 .  A  curved  surface  configuration  (for  the  ground)  is  employed,  which  is  immersed 
in  the  uniform  Cartesian  grid.  A  plane  wave  propagating  through  the  specified  medium  is  shown 
in  Figure  4.6.  Lagrangian  markers  are  also  plotted  for  comparison  purpose.  The  marker  positions 
are  computed  as  given  in  equations  (4.5),  (4.6),  (4.7).  It  can  be  seen  that  the  markers  fail  to 
maintain  the  continuous  surface  when  the  wave  is  diverging,  but  foil  in  the  band  of  contours 
computed  by  WC. 


4.4  Multiple  grids 

One  other  important  use  of  WC  for  the  wave  equation  involves  cases  with  multiple  grids  with 
grid  interfaces.  This  property  plays  an  important  role  when  a  fine  grid  has  to  be  used  around  the 
near  field  (source/antenna)  and  coarser  grid  in  the  far  field.  A  fine  grid  is  typically  used  near  the 
antenna  to  resolve  the  high  field  gradients.  To  propagate  the  wave  over  long  distances,  a  coarse 
grid  scale  must  be  used  in  the  far  field.  The  test  case  used  only  treats  uniform  grids,  and  not  the 
antenna  for  illustration.  If  only  the  discretized  wave  equation  is  used,  with  no  Confinement, 
reflections  are  produced  from  the  grid  interfaces,  unless  special  care  is  taken.  WC  overcomes  this 
problem  and  a  pulse  can  propagate  across  the  interface  with  no  spurious  reflections. 

The  discretized  equation  used  to  solve  the  wave  for  multiple  grids  is 


4.8 


where  ^:=r-r'<  n  = 
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must  be  noted  that  tiie  ratio  —  over  the  entire  domain  is  constant  and  therefore,  the  number  of 
grid  cells  across  the  pulse  remains  constant. 
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•.  ^/=- 


Cl  ^  / 

and  -i-  =  ^  =  const. .  It 


At  the  grid  interface, 


where  values  at  the  interface  are  calculated  as  shown  in  Figure  4.7a.  A  cylindrical  wave  that 
passes  through  the  variable  grid  interface  is  shown  in  Figure  4.7b.  It  can  be  seen  that  there  is 
almost  no  production  of  reflected  waves  from  the  interface.  The  total  amplitude  as  the  wave  is 
passing  through  fine  grid/  coarse  grid  interface  is  shown  in  Figure  4.8.  This  is  checked  for  a  range 
of  grid  ratios  and  presented  in  Figures  4.8a,  4.8b,  4.8c  for  grid  ratios  equal  to  2,  4, 6  respectively. 
Anomalous  numerical  reflections  during  transition  from  the  fine  grid  to  the  coarse  grid  can  be 
suppressed  by  using  variable  confinement  parameters,  which  are  dependent  on  grid  size  and  yet, 

the  ratio  —  can  be  kept  constant. 


5.  CARRIER  FUNCTION 


Developments  that  allow  wave  fronts  to  capture  and  propagate  more  details  are  described  in  this 
chapter.  These  details  are  required  to  compute  interference.  WC  can  be  used  to  generate  constant 
arrival  time  (Eikonal  phase)  surfaces  accurately  by  storing  the  centroid  arrival  time  at  each  grid 
point.  Also,  multiple-arrival  times  are  easily  accomnK)dated,  which  is  complicated  using  Eikonal 
schemes.  Some  recent  developments  in  Eikonal  methods  [2]  can  treat  multiple  arrival  times  but, 
these  methods  require  extra  independent  variables  and  complex  data  management  schemes  are 
used  to  control  memory  requirements.  To  demonstrate  the  ability  of  WC  to  capture  these  surfaces, 
arrival  times  for  the  wave  tonts  generated  from  a  point  source  located  at  (128,128)  on  a  256x256 
grid.  The  parameters  used  for  this  computation  are  v  =  0.23 ,  =  0.2 ,  e  =  0.3 .  The  arrival  time, 

t ,  is  calculated  for  the  passage  of  each  wave  as 


K 


5.1 


It  can  be  seen  in  Figure  5.1a  that  the  contours  of  arrival  times  are  smooth  with  no  stair-casing 
effects,  which  shows  that  there  are  no  plottable  discretization  errors.  The  arrival  times  at  the  grid 
points  in  a  cone  with  a  range  of  angles,  0^29°  -32",  are  compared  to  the  exact  values  for 
validation  purpose.  For  each  i  index,  there  are  several  grid  points  with  varying  j  indices  within  the 
cone  of  ^  =  29" -32".  The  computed  arrival  times  at  each  these  grid  points  is  plotted  against 
exact  values  in  Figure  5.1b.  Assuming  constant  index  of  refraction,  the  exact  arrival  times  are 
computed  as 


=  —  -  5.2 

V 

The  direction  of  propagation,  which  is  the  gradient  of  the  phase  field,  can  also  be  accurately 
computed  as 


5.3 


where  5,  is  the  x-component  of  direction,  and  Sy  is  the  y-component  of  direction.  However,  it  is 

required  that  the  wave  front  has  entirely  passed  the  grid  point  to  accurately  compute  these 
properties.  For  very  small  wavelengths,  the  wave  fronts  and  grid  cells  are  hundreds  of 
wavelengths  wide  and  it  becomes  difficult  to  accommodate  multiple  arrival  times  when  the  waves 
are  close  to  each  other.  To  overcome  that  problem,  the  scalar  field,  ^ ,  is  used  as  a  carrier 
function  that  can  act  as  wave  packet,  which  carries  the  required  details  of  the  propagating 
quantity.  It  is  then  not  necessary  for  the  wave  front  to  entirely  pass  a  grid  point  to  capture  the 
properties.  The  capability  of  WC  is  extended  further  using  this  “carrier’'  approach  and  numerous 
experiments  are  conducted. 


5.1  Propagation  Directions 

The  carrier  approach  is  used  to  accurately  propagate  other  properties  of  the  wave  with  scalar 
function,  ^ .  For  example,  to  propagate  the  directions  in  2D,  the  wave  equation  is  solved  three 

times  for  three  quantities,  ^ ,  ^5,  and  ^Sy  as 


+  ns;  (v  V," )  -  es;  (v^o;  ) 
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where  /  *  1 , 2  and  3  and  is  defined  as, 


^  =4Sy. 

The  updated  directions  are  computed  using  the  normalized  relations. 


s,  = 
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5.6 


In  the  test  computation  described,  the  initial  conditions  for  ^  are  given  as 


with  initial  directions. 


>->0 
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where  the  radius  of  the  wave,  r^j  =  yjii-ioY  +(7“  JoY  origin  at  the  center  of  the  domain 
(%*yo)*  parameters,  v ,  and  e  are  the  same  for  all  three  quantities  and  are  v'  =  0.23, 
//  =  0.2,  and  £  =03.  ^  and  ^2  hehave  as  the  scalar  function  itself  during  wave  collision.  In 
Figure  5.2,  directions  for  the  thin  wave  are  shown  at  4  different  time  steps.  During  interaction, 
directions  temporarily  are  the  mean  values  of  the  waves  that  are  together.  Then,  they  take  their 
original  forms  af^er  a  short  relaxation  time. 


5.2  Edge  Diffractton 

According  to  geometric  optics,  light  travels  in  the  form  of  rays  and  can  bend  in  regions  of  non- 
uniformities  in  the  medium.  In  isotropic  media,  the  rays  are  normal  to  the  primary  wave  fronts 
and  fail  to  accommodate  optical  effects  like  diffraction.  Since,  diffraction  was  not  described  by 
the  classical  ray  theory  of  light  [30],  it  has  been  separately  modeled  using  the  Geometric  Theory 
of  Diffraction.  On  the  other  hand,  Huygens  principle  (wave  optics)  postulates  that  every  point  on 
the  wave  front  act  as  a  point  source  for  further  propagation  and  can  accurately  account  for 
diffraction.  One  such  example  is  diffraction  of  a  plane  wave  from  a  semi-infinite  screen.  The  total 
field  obtained  by  Sommerfeld  [31]  consists  of  incident,  reflected  and  a  cylindrical  diffracted 
wave.  According  to  Huygens  principle,  a  secondary  source  is  created  at  the  edge  and  a  part  of  the 
field  from  that  source  is  transmitted  into  the  non-illuminated  (shadow).  The  incident  field 
interferes  with  the  diffracted  field  and  interference  fringes  are  observed  in  the  illuminated  region. 
At  high  frequency,  the  waves  do  not  bend  and  maintain  sharp  shadows. 

5.2.1  Keller^s  Approximation 

The  diffraction  from  a  knife-edge  can  be  computed  using  the  Geometric  Theory  of  Diffraction 
(GTD)  developed  by  Joseph  Keller  [32]  in  the  1960s.  It  is  an  extension  of  geometric  theory  of 


optics.  This  approximation  will  incorporate  diffraction  effects  into  the  ray  theory  of  light.  Also,  it 
introduces  secondary  sources  that  generate  diffracted  rays  from  the  edges  and  comers  of 
boundaries.  GTD  has  been  successfully  applied  in  acoustics  and  Optics.  Let  us  consider  an 
incident  field  propagating  horizontally  that  encounters  a  knife-edge,  as  shown  in  Figure  5.3b.  The 
fundamental  principle  of  GTD  is  that  ray  propagation  is  a  local  phenomenon.  Furthermore,  all 
fields  of  any  origin  obey  the  laws  of  geometric  optics  locally.  The  same  is  also  true  for  diffracted 
rays.  When  the  incident  field  is  propagated  in  the  direction  normal  to  the  edge,  the  diffracted  field 
is  cylindrical  with  the  edge  as  its  axis.  Now,  the  total  field  at  a  point  is  the  sum  of  all  fields  that 
pass  through  a  point.  An  interference  pattern  is  observed  due  to  superposition  of  the  geometric 
incident  wave  ( )  and  diffracted  edge  wave  ( 


M  =  w'  + 


where 
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5.10 


where  co  is  the  frequency  of  the  wave,  r  is  the  radius  of  the  diffracted  wave  and  D  is  the 
diffraction  coefficient  for  the  knife  edge.  The  diffraction  coefficient,  D.  developed  by  Keller  is. 


D  = 


e 


ir/4 


(2/rky'^  sin  ^ 
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where  0  is  the  angle  between  incident  and  diffracted  rays,  as  shown  in  Figure  5.3b. 


S.2.2  IVC  with  Keller  "s  Approximation 


The  high  frequency  approximation  is  achieved  by  WC  (as  a  pde)  by  restricting  the  dissipation  and 
contraction  to  the  tangential  and  normal  directions  respectively.  The  incident  wave  propagation 
with  a  knife  edge  in  the  path  has  a  sideward  dissipation  when  computed  using  the  equation. 


5.12 


As  explained  in  Section  4.2,  the  nonlinear  term  acts  as  dissipation  in  the  tangential  direction  to 
generate  smooth  wave  fronts.  Due  to  this  effect,  the  arrival  times  start  to  bend  as  the  wave 
encounters  the  edge,  as  shown  in  Figure  5.4.  Even  though  this  looks  like  diffraction,  it  actually  is 
produced  by  dissipation  from  the  confinement  term  to  maintain  the  smooth  field  for  • 


The  scalar  function  can  be  treated  as  a  wave  packet,  which  can  be  used  to  propagate  other 
variables  such  as  direction,  as  demonstrated  in  Section  5.1.  This  allows  WC  to  capture  the 
resulting  incident  wave  produced  by  a  knife  edge  without  numerical  spreading.  The  equation  used 
for  the  high  frequency  approximation  is 


C' = WXj + ] 
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where  s^and  Sy  are  propagated  as  described  in  Section  5.1,  the  components  of  the  normal 
direction  arc  ^y“^x-  "These  waves  obey  the  principles  of  geometric  optics  and  arc 

called  geometric  waves. 

The  arrival  times  of  incident  field,  ,  are  computed  using  the  relation  given  in  equation  (5. 1 ).  It 
can  be  seen  in  Figure  5.5a  that  a  sharp  shadow  boundary  is  generated.  The  edge  wave  is 
propagated  from  a  point  source  placed  on  the  edge  with  the  amplitude  given  by  the  incident  field, 

A* ,  The  arrival  times  of  the  diffracted  wave,  arc  shown  in  Figure  5.5b.  The  total  field  is  now 

calculated  as 


The  interference  pattern  is  computed  using  the  equation  5.1 1  and  compared  to  the  Fresnel  far 

field  approximation  [30].  The  intensity  ratio,  — ,  for  planes,  /=62  and  /  =  68  is  shown  in 

^0 

Figure  5.6a  and  Figure  5.6b  respectively.  It  can  be  seen  that  the  results  agree  closely  with  the 
Fresnel  far  field  approximation. 
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6  CONCLUSION 

The  wave  equation  is  accurately  solved  using  the  newly  developed,  Wave  Confinement  method, 
which  includes  the  realistic  effects  such  as  varying  index  of  refraction,  scattering  from  complex 
objects  and  interference.  Unlike  conventional  Eulerian  schemes,  WC  does  not  suffer  from 
numerical  dissipation  and  the  properties  such  as  total  amplitude  and  centroid  are  conserved.  It 
generates  wave  equation  solutions  as  thin  solitary  waves  that  can  persist  on  the  grid  indefinitely. 
WC  involves  a  nonlinear  term  in  the  propagation  equation,  unlike  the  original  linear  wave 
equation.  However,  this  nonlinear  term  in  the  WC  does  not  produces  amplitude  exchange  or 
phase  shift  during  wave  collisions.  After  collision,  the  waves  take  their  correct  speed  and  form 
after  a  short  relaxation  time.  Also,  the  interactions  behave  in  the  same  way  in  any  dimension. 
This  is  an  important  property  in  capturing  the  effects  of  focusing,  multiple  scattering  and  all  the 
cases  where  many  waves  simultaneously  meet  at  a  point. 

WC  has  been  shown  to  capture  caustic  regions  with  no  additional  steps  in  the  computation 
process.  It  is  also  observed  that  WC  in  combination  with  Ray  Tracing  is  accurate.  This  led  to  a 
new  approximation  called  Eulerian/Lagrangian  Transition.  The  smooth  transition  is  observed 
even  when  there  are  irregular  boundaries.  In  spite  of  having  a  nonlinear  term  in  WC,  the  waves 


arc  observed  to  focus  with  no  shift  in  centroid  positions.  WC  is  also  shown  to  generate  smooth 
scattered  surfaces  during  reflection  from  complex  objects,  preventing  the  need  to  use  complicated 
adaptive  grids.  These  objects  are  immersed  in  a  Cartesian  grid  with  no  grid  refinement  and 
complex  logic.  This  is  due  to  the  confinement  term  being  diffusive  in  the  tangential  direction, 
keeping  the  scattered  field  smooth  with  no  stair-casing  effects.  WC  can  also  accommodate  very 
small  objects  compared  to  the  physical  domain,  and  does  not  require  fine  grids,  unlike 
conventional  Eulerian  schemes.  This  is  an  important  feature  in  electromagnetic  wave  propagation 
in  the  atmosphere  when  complex  surfaces  such  as  terrain,  buildings,  etc.,  have  to  be 
accommodated. 

Another  important  property  of  WC  is  its  capability  to  capture  the  effects  of  varying  index  of 
refraction  efficiently.  This  is  an  important  feature  in  studying  the  effects  of  varying  medium 
properties.  Examples  include  evaporation  or  surface  ducts  (formed  in  atmosphere  due  to 
variations  in  density,  temperature,  etc),  ionospheric  ducts  (formed  in  the  ionosphere  due  to 
variations  in  electron  density).  These  ducts  can  form  caustic  regions,  multiple  scattering  etc.  The 
confinement  term  is  shown  to  treat  such  effects  accurately  and  it  is  validated  using  Lagrangian 
method,  which  is  an  exact  solution  in  high  frequency  approximation.  Wave  propagation  across 
multiple  grid  interfaces  is  also  proved  to  be  accurate  with  no  anomalous  reflections.  This  will 
allow  WC  to  use  fine  grids  around  an  antenna  and  coarser  grids  in  the  far  field.  For  the  tested  grid 
ratios,  the  transfer  of  amplitude  across  the  interface  is  shown  to  be  accurate. 

WC  is  shown  to  accurately  generate  contours  of  constant  arrival  times  as  the  wave  passes  a  grid 
point.  The  scalar  function  propagated  using  WC  can  also  be  used  as  a  carrier  wave  to  propagate 
directions  along  with  the  wave.  This  allows  efficient  computation  of  the  propagation  directions. 
The  properties  propagated  with  the  scalar  function  behave  the  same  way  as  the  scalar  function 
itself  during  multiple  interactions.  Using  the  propagated  directions,  the  wave  solution  in  the  high 
frequency  approximation  can  be  computed.  When  the  wave  encounters  an  edge,  this  will  prevent 
sideward  bending  of  the  waves  into  the  shadow  region  due  to  numerical  diffusion.  Using  Keller's 
GTD  along  with  WC,  the  diffraction  pattern  from  a  knife  edge  is  computed  and  is  shown  to  match 
quite  well  with  the  Fresnel  fiu*  field  approximation. 

The  local  parabolic  method  is  developed,  which  can  compute  propagation  at  any  angle.  It  is  an 
extension  to  die  existing  parabolic  equation,  which  is  confined  to  one  preferred  direction. 
Computations  to  simulate  propagation  at  higher  angles  and  single  slit  diffractions  are  done  and 
validated  using  the  exact  analytical  relations.  Many  computations  are  performed  for  a  number  of 
test  cases  and  validated  to  check  accuracy.  The  numerical  properties  of  WC  have  to  be  well 
understood  before  using  it  for  engineering  applications.  All  the  required  validations  have  been 
performed  for  this  new  method  to  check  the  accuracy  in  actual  complex  conditions.  WC  still 
needs  to  be  validated  for  propagation  and  interference  effects  for  overlapping  waves.  As  the 
pulsewidth  will  be  relatively  large  compared  to  the  wavelength  in  high  fi^uency  propagation, 
work  has  to  be  done  to  compute  such  effects  within  the  wave.  A  more  difficult,  but  widely 
studied  challenge  is  the  inverse  scattering  phenomena.  This  will  be  studied  in  the  future  to  extend 
the  applications  of  WC  further. 


Figure  2.1;  Pulse  simulation  using  equation  (2.30).  The  initial  conditions  for  this 
computation  are  defined  in  equation  (2.29). 


Figure  2.2:  Comparison  to  higher  order  method.  The  pulse  is  propagated  for  100  time 
steps  using  both  WC  and  a  higher  order  method. 
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Figure  3.1 :  Behavior  of  ID  wave  equation  solution 


Figure  3.2:  Computed  pulses  with  and  without  confinement.  The  initial  conditions  are 

given  in  equation  (3.5). 
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Figure  3.3:  Centroid  position.Collison  of  forward  (y  +  vn )  and  backward  {j  -vn 
propagating  pulses  with  each  other  during  propagation 
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Figure  3.4:  Interaction  between  two  pulses  of  different  amplitudes.  The  pulses  undergo 
head-on  collision  during  propagation  for  the  intial  conditions  given  in  equation  (3.10). 


Figure  3.5  :  Diverging  circular  wave  front  solved  using  equation  (3.11) 


Figure  3.6:  Spherical  wave  front  propagation  solved  using  equation  (3.18) 


Figure  4.1:  Inward  propagating  wave  front  computed  using  WC. 
(a)  n  =  0,  (b)  n  =  140,  (c)  n  =  170,  (d)  n  =  250 
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Figure  4.2;  Immersed  boundary  to  treat  complex  boundaries. 


Figure  4.3:  Scattering  from  Infinite  Cylinder 
(a)  n  =  0,  (b)  n  =  200,  (c)  n  =400,  (d)  n  =  600 


Figure  4.4:  Reflection  from  the  boundary  with  shape  of  a  sine  wave 
(a)  n  =  0,  (b)  n  =  200,  (c)  n  =  300,  (d)  n  =  500 
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Figure  4.5:  Varying  index  of  refraction.  The  plane  wave  propagation  across  medium  of 

profile  given  in  equation  4.8 


Figure  4.6:  Plane  wave  propagation  above  curved  surface,  The  medium  has  a  diverging 
profile  given  equation  (4.13)  to  imitate  surface  duct  in  atmosphere 
(a)  n  =0,  (b)  n  =  200,  (c)  n  =400,  (d)  n  =  600,  (e)  n  =  800 


Figure  4.7;  Multiple  grids 
(a)  Computation  of  values  at  interface 
(b)  Circular  wave  front  propagating  across  multiple  grid  interface 
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Figure  4.8:  Amplitude  transfer 

(a)  Grid  ratio  =  2 

(b)  Grid  ratio  =  4 

(c)  Grid  ratio  =  6 
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Figure  5.1:Arrival  times  of  cylindrical  wave  generated  from  point  source 
(a)  Contours  of  arrival  time 
(b)  WC  Vs  Exact  for  0<=29"  -  32" 


Figure  52:  Propagation  of  directions 
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Figure  5.3:  Knife  edge  diffraction 
(a)  Huygens  principle 
(b)  Keller’s  geometric  theory  of  diffraction 


Figure  5.4  :  Sideward  bending  of  arrival  time  contours. 
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Figure  5.5:  Arrival  times  for  knife  edge, 
(a)  Incident  field  /' 

(b)  Diffracted  field 
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